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We study the susceptibilities of quark number, isospin number and electric charge in numerical simulations of 
lattice QCD at high temperature and density. We discuss the equation of state for 2 flavor QCD at non-zero 
temperature and density. Derivatives of InZ with respect to quark chemical potential (fi q ) are calculated up 
to sixth order. From this Taylor series, the susceptibilities are estimated as functions of temperature and p q . 
Moreover, we comment on the hadron resonance gas model, which explains well our simulation results below T c . 



Thermal fluctuations near the critical tem- 
perature (T c ) provide important information for 
the understanding of the properties of the QCD 
phase transition at high temperature and den- 
sity. Baryon number fluctuations are expected 
to diverge at the cndpoint of the first order chiral 
phase transition line. The electric charge fluctu- 
ation in heavy-ion collisions is one of the most 
promising experimental observables to identify 
the critical endpoint. The fluctuations are es- 
timated by the susceptibilities of quark number 
iXq) an d isospin number (xi), which are given 
by Xq — d 2 p/d[ig and xi = d 2 p/d^, where 
^ q = (Hu + Md)/2, Hi = (m« - Md)/2. H u [d] is the 
chemical potential for the u[d] quark. The charge 
fluctuation is xc — Xg/36 + X//4 for \i u = [id- 

The simulation of QCD at non-zero [i q is known 
to be difficult. However, studies based on a Tay- 
lor expansion with respect to /i q at \i q = turned 
out to be an efficient technique to investigate 
the low density regime, interesting for heavy-ion 
collisions. We discussed in pQ the relation be- 
tween the line of constant pressure (energy den- 
sity) and the phase transition line calculating the 



second derivative of pressure. From the study of 
the fourth derivatives [2] , we observed strong en- 
hancement of Xq near T c , suggesting the presence 
of the critical point at finite \i q . In this study, 
we take the calculation of p to 0(/j, q ), and that 
of Xq [xi] to 0(fi q ). This enables us to estimate 
the shift of the peak position of Xq m t ne C^i A*g) 
plane. Also, by calculating the ratio of the Tay- 
lor expansion coefficients, the application range of 
the hadron resonance gas model |3| is discussed. 

Quark gluon gas, hadron resonance gas 

We define the expansion coefficients c„ and 
by p/T 4 = (lnZ)/(UT 3 ) ee En=o^ q /T) n and 
Xq[xi\/T 2 = En= 2 n(n - 1 )c„[4](^/T)- 2 for 
= Hd- We expect the equation of state to 
approach that of a free quark-gluon gas (Stefan- 
Boltzmann (SB) gas) in the high temperature 
limit. The coefficients in the SB limit for N{ = 2, 
pi = are well known as C2 = c\ = 1 , C4 = c\ = 



1/(2tt 2 ) and c„ 



for n > 6. 
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On the other hand, in the low temperature 
phase QCD is well modelled by a hadron res- 
onance gas. If the interaction between these 
hadrons can be neglected, the pressure is obtained 
by summing over the contributions from all res- 
onance states of hadrons. The contribution to 
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Figure 1. The ratio of Taylor expansion coeffi- 
cients for Xq (upper) and \i (lower). 

p/T 4 from a particle which has mass rrii, baryon 
number Bi and third component of isospin num- 
ber I^i is given by 
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where z = exp[(3£?i/i 9 + 2/3^/) /T], is the 
modified Bessel function, and 77 is 1 for mesons 
and —1 for baryons. Moreover, since rrii/T 3> 1 
for all baryons and K2(x) « y / ir/2x exp(— x), the 
first term of I is dominant in the baryon sector. 
Therefore, the pressure can be written as p/T 4 — 
G(T)+F(T) cosh(3/i,/T) for m = 0, where F(T) 
is the baryonic component of p/T at fi q — and 
G(T) = G sing {T) + Gtrip(T) is the mesonic part 
which has isosinglct and isotriplet components 3 . 
Similarly, we obtain Xq /T 2 = 9F{T) cosh(3fi q /T) 
and xi/T 2 = G 7 (T) +F / (T)cosh(3^ 9 /T). Here, 
the mesonic component for Xq 1S zero because 
mesons (Bi = 0) are independent of \x q . There- 
fore, C4/C2 = 3/4, c 6 /c 4 = 3/10 and c\/c\ = 3/10 
in the region where the non-interacting hadron 
resonance gas provides a good approximation. 
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Figure 2. Pressure at /J, q = as a function of T. 

We investigate these coefficients. Simulations 
are performed on a 16 3 x 4 lattice for 2 flavor 
QCD. The p4 improved action [2| is employed 
at ma = 0.1, which gives mps/my — 0.7. The 
number of configurations is 1000-5000 for each T. 
We use the random noise method 12. The re- 
sults for c„+2/cn and c^ +2 /c^ are shown in Fig.^ 
We find that these results are consistent with the 
prediction from the hadron resonance gas model 
for T/T c < 0.96 and approach the SB values, i.e. 
c 4 /c 2 = C4/C2 = 1/2-7T 2 , c 6 /c 4 = c\/c\ = 0, in 
the high temperature limit. These results sug- 
gest that the models of free quark-gluon gas and 
hadron gas seem to explain the behavior of ther- 
modynamical quantities well except in the narrow 
regime near T c . 

Mesonic and baryonic contributions 

Moreover, if the data is consistent with the 
hadron resonance gas model, we can separate xi 
into the contributions from isotriplet meson and 
baryons, since G 1 + F 1 = 2c/ 2 and 9F 1 /2 = 12c{. 
Isosinglet mesons do not contribute to xi because 
Isi = 0. Also, 9F — 2c2 from Xq- Furthermore, 
for the case of ^S> T, the first term (1 = 1) 
in Eq.Q dominates; this approximation holds 
for our present simulations where m^/T ~ 4-5 
but will no longer be valid with physical quark 
mass values. In this approximation, the ratio of 
d 2 {p ml /T 4 )/d(^/T) 2 and Pml /T 4 equals (2I 3l ) 2 
for I 3i = {—1,0, 1}. We thus find for the contri- 
bution of the triplet meson part, G tr ip = 3G 7 /8. 
We then can estimate the contribution of the 
isosinglet mesons to p/T , G s j ng , by comparing 



3 



F + Gtrip with the total pressure. In Fig. |3 we 
plot the baryon component F = 2c 2 /9 (square), 
the triplet meson component G tr i p = ic\l^ — c| 
(circle) and the total (diamond) . The solid line is 
the pressure obtained with the integral method in 
[5] . The dashed lines in Fig. [2] represent the pre- 
dictions from the hadron resonance gas model for 
the mass parameter of our simulations. The res- 
onance states for this quark mass are adjusted by 
the method described in The simulation re- 
sults for p/T 4 obtained with the integral method, 
together with the triplet meson and baryon con- 
tributions, are surprisingly well reproduced by 
this model calculation. The result for the singlet 
meson contribution may explains the difference in 
pressure obtained from the integral method and 
that from the Taylor expansion which neglects the 
singlet part. 

Susceptibilities at /i q ^ 

Next, we calculate quark number and isovec- 
tor susceptibilities, using three methods, in a 
range of < Hq/T < 1. The data connected 
by solid lines in Fig.|31are obtained by Xg/T 2 = 
2c 2 + Uaifiq/T) 2 + 30c 6 (fi q /T) 4 and the corre- 
sponding equation for \i- Dot-dashed lines are 
from the hadron resonance gas model, using F, F 1 
and G 1 . This turned out to be a good approxima- 
tion for T/T c < 0.96. Dashed lines are calculated 
by the reweighting method with an approxima- 
tion in pp. Here the truncation error is 0(/ig) 
but the effect from higher order terms of \x q is par- 
tially included. The sign problem in the reweight- 
ing method is serious when the complex phase 
fluctuations of the quark determinant are large 
at large n q /T. We omitted data for which the 
standard deviation of the complex phase is larger 
than 7r/2. The difference among the three results 
is caused by the approximation in the higher or- 
der terms of [i q jT. 

Since the statistical error of eg is still large 
near T c , the peak of \q seems to lose its sta- 
tistical significance. However, as seen in Fig. ^ 
cq changes its sign at T c . This means the peak 
position of \ q moves left, which is correspond- 
ing to the change of T c as a function of [i q . 
T c (p q /T = 1)/T c (fj, q /T = 0) in [T] is about 0.93. 
Moreover, Xq increases with higher orders of the 




Figure 3. Quark number (upper) and isovector 
(lower) susceptibilities as a function of T for each 
fixed Hq/T. To is T c at /i g = 0. 

expansion for T <T C which confirms the existence 
of a peak. This suggests the presence of a criti- 
cal endpoint in the (T, fj, q ) plane. On the other 
hand, \i m Fig.[3]does not show any singular be- 
havior. This is consistent with the sigma model 
prediction that only isosinglet degrees of freedom 
become massless at the critical endpoint. 
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